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Power transmission in one-dimensional nonlinear magnetic metamaterials driven at one end is 
investigated numerically and analytically in a wide frequency range. The nonlinear magnetic meta- 
materials are composed of varactor-loaded split-ring resonators which are coupled magnetically 
through their mutual inductances, forming thus a magnetoiductive transmisison line. In the lin- 
ear limit, significant power transmission along the array only appears for frequencies inside the 
linear magnetoinductive wave band. We present analytical, closed form solutions for the magne- 
toinductive waves transmitting the power in this regime, and their discrete frequency dispesion. 
When nonlinearity is important, more frequency bands with singificant power transmission along 
the array may appear. In the equivalent circuit picture, the nonlinear magnetoiductive transmis- 
sion line driven at one end by a relatively weak electromotive force, can be modeled by coupled 
resistive-inductive-capacitive (RLC) circuits with voltage-dependent capacitance. Extended numer- 
ical simulations reveal that power transmisison along the array is also possible in other than the 
linear frequency bands, which are located close to the nonlinear resonances of a single nonlinear RLC 
circuit. Moreover, the effectiveness of power transmission for driving frequencies in the nonlinear 
bands is comparable to that in the linear band. Power transmission in the nonlinear bands occurs 
through the linear modes of the system, and it is closely related to the instability of a mode that is 
localized at the driven site. 
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Keywords: nonlinear split-ring resonators, magnetic metamaterials, driven linear waves, nonlinear power 
transmission 



I. INTRODUCTION 

The metamaterial concept usually refers to a periodic 
arrangement of artificially structured elements, designed 
to achieve advantageous and/or unusual electromagnetic 
properties compared to those of naturally occuring mate- 
rials i42i] . For example, the first metamaterial ever real- 
ized [40| , exhibited negative refraction index in a narrow 
frequency band around 5 Gigahertz. The negative re- 
fraction index property requires that both the dielectric 
permittivity and the magnetic permeability are simul- 
taneously negative. The most widely used elements for 
constructing a negative refractive index metamaterial are 
the electrically small resonant 'particles' called split-ring 
resonators (SRRs) and metallic wires. While the array of 
wires gives a negative permittivity below the plasma fre- 
quency, the array of the SRRs gives a negative magnetic 
permeability above their resonance frequency. Thus, the 
refraction index can be negative in a narrow frequency 
band. Since their discovery, metamaterials have been the 
target of intensive study. Although the developement of 
metamaterials at microwave frequencies has progressed 
to the point where scientists and engineers are now pursu- 
ing applications, research on metamaterials that operate 
at higher frequencies is still in an early stage [1^, [la, . 

There is a wide subclass of metamaterials, that exhibit 
significant magnetic properties and even negative mag- 
netic permeability at Terahertz and optical frequencies 
[13l |25| . l5l| , even though they are made of non-magnetic 



materials. Magnetism at such high frequencies is partic- 
ularly important for the implementation of devices such 
as tunable lenses, adaptive mirrors, etc., since there are 
only a few natural materials that respond magnetically 
above microwaves. Moreover, magnetism in those natu- 
ral materials is usually weak and within narrow frequency 
bands, limiting thus their use in possible Terahertz de- 
vices. The most common realization of such a magnetic 
metamaterial is composed of periodically arranged SRRs 
in a one- or two-dimentional lattices. In its sirnplest ver- 
sion an SRR is just a highly conducting [sl, [H, US] or 
superconducting Q ring with a slit. The SRR structure 
was originally proposed by Pendry [30| , who theoretically 
predicted that a periodic SRR-based structure could ex- 
hibit a negative magnetic permeability frequency band 
when irradiated by an alternating electromagnetic field 
of appropriate polarization. Pendry also realized that the 
SRR structure has considerable potential for enchancing 
nonlinear phenomena. When the SRR dimensions are 
much smaller than the wavelength of the incident field, 
the SRR can be regarded as an electric circuit consisting 
of inductive, capacitive and ohmic elements connected in 
series. The ring forms the inductance L, while the slit 
can be considered as a capacitor of capacitance C. The 
ohmic element R may represent all possible losses of the 
structure. 

The next step was naturally to built nonlinearity in the 
SRR structure. Nonlinearity implies real-time tunability 
and multistability that is certainly a desired property of 
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possible future devices. The SRRs can become nonlinear 
either by the insertion of a strongly nonlinear dielectric 
or a nonlinear electronic component (i.e., a varactor) 
3^I13,[1q| into their slits. Both ways would lead to effec- 
tively field dependent magnetic permeability for an array 
of nonlinear SRRs. Recently, the dynamic tunability by 
an external field for a two-dimensional array of varactor- 
loaded SRRs was demonstrated experimentally [s^l ■ The 
SRRs in such an array are coupled through magnetic 
and/or electric dipole-dipole interactions, whose strength 
depends on the relative orientation of the SRRs and their 
slits in the array [1, [13, HH, SI] . 

In the equivalent circuit picture, the dynamics of a 
periodic nonlinear SRR array irradiated by an electro- 
magnetic field, can be described by a set of coupled and 
driven ordinary differential equations with on-site non- 
linearity. The discreteness, which is inherent in those 
SRR-based magnetic metamaterials, along with the non- 
linearity and the weak coupling between their elements, 
allow for the generation of nonlinear excitations in the 
form of intrinsic localized modes or discrete breathers 
(4]| . These nonlinear modes appear generically in dis- 
crete and nonlinear extended systems. Recent theoreti- 
cal work in one- and two-dimensional nonlinear SRR lat- 
tices have demonstrated the existence and the stability of 
discrete breathers both in energy-conserving and dissipa- 
tive systems [1, 0, [l^ [l^ . It has also demonstrated that 
breathers may be formed spontaneously through modu- 
lational instability in binary nonlinear magnetic meta- 
materials [1^ [12, [i^ . Moreover, domain- wall excitations 
[35I I and envelope solitons [H, [2, [3 El] may as well be 
excited in those systems, which seem to be stable even 
in the presence of noise [49] . 

The SRR-based magnetic metamaterials support a 
new kind of electromagnetic waves, the magnetoinduc- 
tive waves, which exhibit phonon-like dispersion curves, 
and they can transfer energy along the array [s^, 13 E3| ■ 
It is thus possible to fabricate a contact-free data and 
power transfer device, a magnetoinductive transmission 
line, which make use of the unique properties of the 
magnetic metamaterial structure, and may function as 
a frequency-selective communication channel for devices 
via their magneto-inductive wave modes j44| . In the 
present work we investigate the power transmission along 
a one-dimensional varactor-loaded SRR-based magnetic 
metamaterial, in the two possible geometries, which is 
driven at one end by a sinusoidal power sourse. In the 
linear limit, energy tranfer along the array occurs only for 
driving frequencies in the linear magnetoinductive wave 
band. However, the nonlinearity could generate more fre- 
quency bands where efficient power transmission along 
the array is possible. Power transmission in chains of 
coupled anharmonic oscillators for driving frequencies in 
the band gap of the linear spectrum has been recently 
investigated, and that effect is referred to as self-induced 
transparency (28j or supratransmission (3]. For frequen- 
cies inside the linear magnetoinductive wave band the 
power is transmitted along the array from the magne- 



toinductive wave modes, slightly modified by the pres- 
ence of nonlinearity. However, for frequencies in the pass- 
bands resulting from nonlinearity, the power is transmit- 
ted along the array from nonlinear, breather-like exci- 
tations. In the next Section we derive the model equa- 
tions for a varactor-loaded SRR-based magnetic meta- 
material, where the nonlinearity has the form a polyno- 
mial expansion on the main variable, like that investi- 
gated in ■ We also present bifurcation diagrams for 
the shortest possible end-driven array, i.e., an array with 
only two varactor-loaded SRRs, from which only one is 
directly driven. In Section HI we obtain analytical solu- 
tions for the linear magnetoinductive waves of the end- 
driven, finite SRR array when losses are neglected, along 
with their discrete dispersion relation. In Section IV we 
present numerical results for power transmission and dis- 
cuss their dependence on the model parameters for a 
short, end-driven varactor-loaded SRR array with and 
without an absorbing boundary at the non-driven end. 
We finish with concluding remarks in Section V. 



II. VARACTOR-LOADED SRR ARRAY MODEL 

Consider a ring-shaped or a squared split-ring res- 
onator with a hyperabrupt tuning varactor mounted 
onto its slit. The varactor type could be selected 
from a large variety of available varactors, according 
to the needs of the experiment. In recent experimen- 
tal works on varactor-loaded SRR-based metamateri- 
als [nl, m, m, 113' the varactor selected was a Sky- 
works SMV1231-079, whose voltage-dependent capaci- 
tance C{Ud) {Ud is the voltage across the diode) is given 
by 



CiUo) = Co 1 



Ud 



-M 



(1) 



where Co is the DC rest capacitance. Up is the intrinsic 
potential, and M is a parameter. The values of those pa- 
rameters are provided by the manufacturer SPICE model 
for that varactor to be 



M = 0.8, Co = 2.2 pF, Up = 1.5 V. 



(2) 



From Eq. ([T]) we can determine the voltage-dependence 
of the normalized charge 



Qi 



CoUp ' 



(3) 



with Qu being the charge in the diode, as a function of 
the normalized diode voltage u — Ud /Up as 



1 



q = 



-M + 1 



l-il-u) 



-A/+1 



(4) 



Assuming that Ud < Up, the earlier equation can be 
solved for the diode voltage as a function of the charge, 
i.e., u = u{q) 



1- [1 -q{-M + 1)] (-^^+1) 



(5) 
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When |a;| = M+l)| << 1, the voltage w in the earlier 
equation can be expanded in a Taylor series inpowers of 
q around zero. Thus, neglecting term of order 0{q'^) or 
higher, we get 



aq 



(6) 



where, by using Eq. ([2]), 
M 



(3 = -M{2M - 1) = +0.08. (7) 
6 



The expression Eq. © can then be used in the volt- 
age equation, obtained from Kirchhoff's voltage law, 
which represents the varactor-loaded SRR as an effective 
resistive-inductive-capacitive (RLC) circuit with voltage- 
dependent capacitance 



cPQl 



dt 



(8) 



where 



UpU = Ud = Up 



CqUp 



Q 



D 



CalJp 



CoUp 



(9) 

, £{t) is the induced electromotive force of amplitude So 
and frequency w resulting from an applied electromag- 
netic field with the same frequency, of the form 



£{t) — £q sin(a;t). 



(10) 



where L is the inductance of the ring, R is the sum of 
the Ohmic resistance of the SRR ring and the series re- 
sistance of the diode, and t the temporal variable. The 
eigenfrequency wq of the circuit described by Eq. ([8]) is 
given, in the linear limit without losses and driving, by 
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Wo 



(11) 



That approximate model of the varactor-loaded SRR is 
found to be adequate for varactor voltages not exceeding 
0.5 V {u < 1/3 in normalized units), when the dissipative 
current of the varactor can be neglected [s^] • 
Using Eq. ([3]) and the relations 



Wo 



Up' 



(12) 



Eq. (IS]) can be written in the normalized form 

, , 2 , « 3 , '^(1 
-j^ + q + aq + (jq +7^ 



: £0 sin(f2r). 



where 
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= R^/Co/L 



(13) 



(14) 



is the normalized loss coefficient. The varactor-loaded 
SRR is a nonlinear oscillator, which exhibits multista- 
bility, hysteretic effects, and secondary resonances, when 




FIG. 1: The current amplitude io as a function of the normal- 
ized frequency for a varactor-loaded SRR modeled by Eq. 
(|13p . Red curve: frequency increases; black curve: frequency 
decreases. 



driven by a sufficiently strong external field. A typical 
current amplitude io - frequency Q curve for the equiva- 
lent circuit model Eq. (fT3|) . like that shown in Fig. 1, ex- 
hibits all those characteristics. Notice the hysteresis loop 
at normalized frequency around 1.1, close to the reso- 
nance frequency of the linear system, and the very strong 
resonance at its second harmonic, i.e., at ft ~ 2, where 
a smaller hysteresis loop also appears. Furthermore, we 
observe a subharmonic resonance at $7 ~ 0.48, as well as 
weaker resonances at f2 ~ 0.63, 0.33 and fl ~ 0.25. For 
frequency intervals lying between the boundaries of the 
hysteresis loops there are two different states of the os- 
cillator, with low and high current amplitude, which are 
simultaneously stable. 
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FIG. 2: Schematic of an SRR array in the planar geometry 
(upper) and in the axial geometry (lower). 



Consider a one-dimensional periodic array of identical 
varactor-loaded SRRs, which are coupled magnetically to 
their first neighbors through their mutual inductances. 
The array can be formed in two possible geometries, ac- 
cording to the relative orientation of its elementary units 
(i.e., the SRRs) in the array. As can be seen in Fig. 2, 
the SRRs can be arranged either in the planar geome- 
try, where all of them are lying on the same plane, or in 
the axial geometry, where the axes of all the SRRs are 
lying on the same line. The orientation of the incident 
electromagnetic field, which is also shown in the figure, is 
such that the magnetic component can excite an electro- 
motive force in each SRR. Below we consider arrays that 
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are driven only at one end (let us say the left end). Thus 
only one SRR receives energy directly from the driver. 
However, the energy can be transmitted along the array 
due to the interaction between its elements. Including 
the coupling of each SRR to its nearest neighbors, we get 
for the charge Qo^n of the n— th varactor the dynamic 
equations 



L- 



dln 

dt 



RIn 



M 



din-l 

dt 



M- 



.dl 



n+l 



dt 



where n = 1,...,A^, with N being the total number of 
SRRs in the array, — dQu^n/dt is the current flow- 
ing in the n— th SRR, UpUn is the approximate voltage 
across the n— th varactor, and M (not to be confused 
with the varactor parameter) is the mutual inductance 
between two neighboring SRRs. Using the relations from 
Eq. the coupled nonlinear Eqs. can be written 

in normalized form as 

d^ 

— (A qn-i + 9„ + A Qn+i) + qn + aql + PqI = 

+£osin(17T)(5„,i, (16) 

where A = M/L is the coupling coefficient between neigh- 
boring SRRs. The function 5„ i is unity for n = 1, while 
it is zero for any other n. The coupling parameter A may 
assume both positive and negative values, corresponding 
to planar and axial geometry of the varactor-loaded SRR 
array, respectively. 

For a particular array, the values of the coupling and 
loss coefficients A and 7, respectively, can be estimated 
from its geometrical and material parameters. Specifi- 
cally, the mutual inductance M, the ring inductance L, 
and the Ohmic resistance Rsrb. of the ring can be cal- 
culated from the following expressions 



M 
L 

RSRR 



/io- 



41)3 ' 



Mo a 
2a 

ahS ' 



In 



16a 



1.75 



(17) 
(18) 
(19) 



where /io is the magnetic permeability in vacuum, a is 
the average radius of each SRR, h the diameter of the 
wire of each SRR, D is the center to center distance be- 
tween neighboring SRRs, and a and S is the conductivity 
and skin depth of the SRR wire. Eqs. (II71)-([I11) refer to 
circular SRRs with circular cross-section. However, they 
can also be used for circular SRRs with square cross- 
section by defining an equivalent diameter of the SRR 
wire h = y^Adw/n, with d and w being the depth and 
width, respectively, of the SRR wire. For a ~ 3.7 mm, 
Z? ~ 11 mm, h ~ 0.5 mm, we get from Eqs. (jl7p and 
(HH), respectively, that M ~ 0.14 nH and L ~ U nH, 



and consequently A ~ 0.01 . For 5 ~ 2.2 x lO^^ mm 
and a ~ 0.5 x 10^ (ilcm)"^ (for an SRR made of copper 
wire at ^ 1 GHz at room temperature) we get from Eq. 




FIG. 3: (a) Energy Etot bifurcation diagram for varying fre- 
quency Q. for an asymetrically driven varactor-loaded SRR 
dimer with 7 = 0.001, eo = 0.3, and A = -1-0.01. (b) Energy 
fractions /i and /2 of the first (driven) and the second SRR, 
respectively, of the dimer in the low energy state, (c) Energy 
fractions /i and /2 of the first (driven) and the second SRR, 
respectively, of the dimer in the high energy state. 



p9)) that RsRR — 0.13 Q.. The diode series resistance Rs 
of the varactor used in Ref. iSO] is around ~ 2 Jl; that 
is, most of the Ohmic resistance of the varactor-loaded 
SRR, R = RsRR + Rs = 0.13-^2.00 = 2.13 n comes from 
the varactor. By substituting i? ~ 2.13 f2 in the second 
of Eq. (HH), along with L = U nH and Cq = 2.2 pF, we 
get that 7 = 0.027. The resonance frequency of a single 
varactor-loaded SRR is fo ~ {27r^/L(J^)-^ ~ 0.9 GHz. 
In the numerical calculations below we used the esti- 
mated value of A that indicates weakly coupled SRRs. 
However, we use a smaller value of 7, which could be 
achieved in practice by using a metalic wire with higher 
conductivity and an ultra-low resistance varactor with 
similar capacitance- voltage relation. 

The shortest end-driven varactor-loaded SRR array is 
that comprised of two elements. Since only one of them is 
driven by the external field, we refer to this system as the 
asymmetrically driven nonlinear dimer. The increase of 
degrees of freedom leads to interesting and complicated 
behavior, that is marked by the appearance of chaos. In 
Fig. 3, the bifurcation diagram of the total energy Etot of 
the dimer with varying normalized frequency D, is shown 
for positive coupling A. For negative coupling A with 
equal strength we get very similar results that share a 
number of characteristic features. First of all, multista- 
bility is again present, for normalized frequencies in an 
interval around 0.9. There, two simultaneously sta- 
ble energy states, with high and low energy, co-exist, as it 
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is shown in Fig. 3a. The energy difference between those 
states increases condiderably with increasing frequency. 
Moreover, low energy state becomes chaotic for frequen- 
cies around ~ 1.9. It is also interesting to see how the 
total energy is divided between the two elements. The 
energy fractions /i = Ei/Etot and /2 = E2/Etot, with 
El and E2 being the energies of the first (driven) and 
second SRR of the dimer, are shown both for the high 
and the low energy states in Fig. 3b and 3c as a func- 
tion of the normalized frequency fi. When the dimer is 
in the low energy state, most of the energy is concen- 
trated in driven SRR whose energy fraction is close to 
unity (correspondingly, the energy fraction of the other 
SRR is close to zero). However, for frequencies around 
the resonances, energy can be transfered easily from the 
one SRR to the other, and then the energy fractions of 
the SRRs in the dimer attain comparable values. In the 
chaotic regions, the energy is tranfered irregularly from 
one SRR to the other and the energy fractions span the 
whole interval between zero and one. In the high energy 
state (Fig. 3c), the energy fractions are almost equal for 
frequencies above the linear resonance frequency, n 1. 
Below that frequency, the energy is either concentrated in 
the driven SRR (for frequencies far from resonances), or 
the values of the energy fractions are comparable (close 
to resonances. 



III. DRIVEN LINEAR SOLUTIONS 

For a weakly driven array, the nonlinear terms which 
are proportional to a and /3 are not so important and 
they can be neglected. Then, the system of Eqs. (fTT 
reduces to a linear one 



the structure, i.e.. 



{Xqn^i - qn + Ag„+i) -I- q„ 
= -1—, l-£o sm(i2r)d„,i. 

(XT 



(20) 



The dispersion relation of the linear system can be ob- 
tained by the substitution q„ = A cos(Avn — fir) in the 
absence of losses and applied field, i.e., for eg = and 
7 = 0, which gives the linear dispersion of magnetoiduc- 
tive waves 



1 



^1 2 A cos(k) 



(21) 



where k is the (normalized) wavevector (— tt < k < tt). 
Eq. (PT|) defines a frequency band of width Afi ~ 2A 
that is bounded by a minimum and a maximum frequency 
^min — l/Vl + 2A and ^max — 1 / V 1 — 2A, respectively. 

However, it is also possible to calculate the dispersion, 
expressed as a series of resonant frequencies, and the ex- 
act form of the linear magnetoinductive modes for the 
finite system in the absence of losses (7 = 0). For finite 
systems, Eq. (|20|) should be implemented with free-end 
boundary conditions, to account for the termination of 



90 = = 0. 



(22) 



By substituting qn — Qn sin(5lT) into Eq. ([20| . the sta- 
tionary equations can be written as 



sQn~l + Qn + sQn +1 = K, 



(23) 



where 



£0 



jSn.i, (24) 



1 - ' 1 _ f]2 

or, in matrix form 

Q = kS-1Ei, (25) 

where Q and Ei are iV— dimensional vectors with com- 
ponets Qi and Si^i, respectively, and is the inverse 
of the N X N coupling matrix S. The latter is a real, 
symmetric tridiagonal matrix that has diagonal elements 
equal to unity, while all the other non-zero elements are 
equal to s. The need to find the inverse of tridiagonal 
matrices like S arises in many scientific and eng ineering 
applications. Recently, Huang and McCoU [l^] related 
the inversion of a general trigiagonal matrix to second 
order linear recurrences, and they provided a set of very 
simple analytical formulae for the elements of the inverse 
matrix. Those formulae lead immediately to closed forms 
for certain trigiagonal matrices, like S [2l|. The compo- 
nents of the Q vector can be written as 

= /^(s-i)^^, (26) 

where (^^^^^ is the (zl)— element of the inverse cou- 
pling matrix S~^, whose explicit form is given in f2l|. 
Then, the solution of the linear system Eq. (|20p and 
(j22p for any driving frequency and finite N is given by 



sm[{N - n + 1)0'] ^. 
sin[(A^ + 



Y)^Bin(17r) 



' = cos-^ -3 , (27) 



2|s 



for s > +1/2 and s < -1/2, and 



smh[{N-n + l)e] . 



^lui±43M!, (28) 



2 s 



for -1/2 < s < +1/2, where 



1 f \s\ 



(29) 



Note that the QiS are uniquely determined by the pa- 
rameters of the system. 



n 




FIG. 4: Density plots of the Qn's in a stationary state, on the 
site number n - normalized frequency plane calculated from 
the analytical expressions Eqs. (|27p and (|28p for A — +0.01 
and £o ~ 0.3. The vertical white bands indicate the positions 
of the resonances. 



From the analytical solution in the linear magentoin- 
ductive wave band, Eq. (I27p . which corresponds to ei- 
ther s > +1/2 or s < —1/2, we can infer the resonant 
frequencies from zeroing the denominator, i.e., by setting 
sin[(iV + 1)9'] — 0. Then we see that the solution has a 
resonance for 



1 



2 cos 



(N + l) 



(30) 



where m is an integer (m = 1, N). The resonant fre- 
qquencies are obtained by solving the first of Eq. (p4)) 
with respect to il, and substituting the resonant values of 
s = Sm from Eq. (j30]) . Then we get the discrete dynamic 
dispersion relation of linear magnetoinductive waves in 
the varactor-loaded SRR array, the discrete analogue of 
Eq. (121]), as 



1 



l-2Acos(^ 



(31) 



where m is the mode number (to — 1, N). In Fig. 4, a 
density plot of the Q^'s on the n — fl plane, calculated an- 
alytically from Eqs. p7)) and p8|. is shown for positive 
coupling coefficient A. In the frequency interval shown, 
the QnS have rather large values because, due to the nar- 
row band, all frequencies are very close to a resonance. 
Moreover, the analytical solutions presented above have 
been obtained in the absence of losses, so that there is 
nothing to prevent the Q„'s from increasing indefinitely 
close to a resonance (and to go to infinity exactly at reso- 
nance). The white vertical regions in these figures corre- 
spond to frequency intervals around a resonance, where 
the QnS have very large values. Thus, those regions in- 
dicate the frequencies of the resonances. We can also 
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FIG. 5: Typical profiles calculated analytically from Eqs. 
1271) and 1211 for A = +0.01, eo = 0.3, and five diflferent 
frequencies: (a) Q. = 0.989; (b) Q. = 0.995;(c) Q = 1.000; (d) 
Q — 1.005; (e) = 1.011. The frequencies in (a) and (e) lie 
outside the linear magnetoinductive wave band. 



see in Fig. 4 that the solutions close to the lower bound 
(S^min — 0.99) of the linear magnetoinductive wave band 
exhibit very high amplitudes. For negative coupling co- 
efficient with equal strength the high amplitude solutions 
appear close to the upper bound {^max — 1.01) of the 
linear magnetoinductive wave band. 

In Fig. 5 the profiles of the analytical solutions for five 
different frequencies, both inside and outside the linear 
magnetoinductive wave band are shown for positive A. 
The frequencies in Figs. 5a and 5e lie outside the lin- 
ear band, so that the Qn's decrease exponentially with 
increasing site number n. 



IV. POWER TRANSMISSION ALONG THE 
ARRAY. 

In the following we focus on the power transmission 
along a varactor-loaded SRR array with = 10 ele- 
ments, which is driven at the left end (n = 1). The 
average power dissipated in the n— th SRR is defined as 



Pn^R< II{t) >, 



(32) 



where <>ro denotes time-average over tq time-units. 
The power density P„ can be normalized to Pq = IqR, 
where Iq = ujq Cq Up and R the value of resistance that 
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results from the value of the dissipation coefficient used 
in the simulations. For 7 — 0.001 we have that i? ~ 0.06 
and Pq — 28 fiW. It is convenient to express P„ in dBm, 
according to the relation 



PnidBm) = 20 log {<ii> 



Poiin W) 
1 mW 



4.341n(0.028 < il >, 



(33) 



where in is the normalized current in the n— th SRR. 
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FIG. 6: Average power density plot (in dBm) on the n — f 
plane, for A = +0.01, 7 = 0.001, and £0 = 0.3, in the linear 
regime (a = 0, P — 0). 



The system of Eqs. is integrated with a 4th 

order Runge-Kutta algorithm with fixed time-stepping, 
and boundary conditions qo = Qn+i = to account for 
the termination of the structure. The initial conditions 
were set to zero. The time allowed for the elimination 
of transients in each value of the frequency is 2000 Tq ( 
To = 27r/a;o), while the time interval tq over which the 
average power is calculated is typically 1000 Tq. In order 
to convert time intervals to physical units one has to mul- 
tiply by the inverse of the resonance frequency Jq^, which 
gives To — 2.2 /zs. In the linear regime, i.e., a = and 
/3 = 0, the average power P{dBm) density plot on the site 
number n - frequency / plane shown in Fig. 6 exhibits 
the expected behavior; that appreciable power transmis- 
sion only occurs for frequencies in the linear band, i.e., 
around / = 0.9 GHz. However, the activation of the 
nonlinear terms reveals that even a relatively low driv- 
ing amplitude eo = 0.3 is sufhcient for nonlinear power 
transmission bands to bew formed. Thus, power can be 
transmitted not only for frequencies in the linear band, 
but also in the otherwise forbidden frequency regions. As 
it can be observed in Fig. 7, a number of transmission 
bands have appeared due to the nonlinearity. The ar- 
ray thus becomes transparent in power transmission at 
frequency intervals around the nonlinear resonances of a 
single varactor-loaded SRR oscillator. That self-induced 
transparency is a robust effect and survives with consid- 
erably increased losses, as it is shown in Fig. 8, where 
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f (GHz) 



FIG. 7; Average power density plot (in dBm) on the n — / 
plane, for A = +0.01, 7 = 0.001, a = -0.4, /3 = 0.08, and 
£0 = 0.3. 



the loss coefhcient has been increases by an order of mag- 
nitude (from 0.001 to 0.01. In the density plots shown 
above, significant power transmission is indicated by red 
color. In both Figs. 7 and 8, nonlinear power trans- 
mission occurs at / ~ 0.3 GHz, f ^ 0.45 GHz, and 
f ^ 1.8 GHz. The latter nonlinear band, around the 
second harmonic of the resonance frequency of the linear 
system, seems to be the wider of all. 
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FIG. 8: Average power density plot (in dBm) on the n — f 
plane, for A = +0.01, 7 = 0.01, a = -0.4, /3 = 0.08, and 
£0 = 0.3. 

When losses are small, the appearance of more non- 
linear bands is possible, like the one seen just above the 
linear band in Fig. 7, at / ~ 1 GHz. That transmission 
band is due to an interplay between nonlinearity and ge- 
ometrical (lattice) resonances, and it is very sensitive to 
changes in the parameters of the model. Note that in the 
density plots presented above the frequency is given in 
natural units, using the relation f = f^Vl = 0.957 GHz. 
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It is very illuminating to see the time evolution of the 
current flowing in each SRR of the array for some fre- 
quencies in the transmission bands. We have chosen two 
different frequencies, one in the linear band and one the 
high frequency band which is around the second har- 
monic of the linear resonance frequency. The time evo- 
lution of the currents in the first, fourth, seventh, and 
tenth (last) SRR are shown for / = 0.888 GHz and 
/ = 1.78 GHz in Figs. 9a and 9b, respectively. The 
currents oscillate periodically in all SRRs, although with 
different amplitudes and relative phases with respect to 
the driver. This figure suggests that power transmis- 
sion for frequencies in the forbidden band gap occurs 
through the linear modes of the system. For frequen- 
cies outside the linear and nonlinear bands, the currents 
exhibit similar oscillations, although their amplitude be- 
comes vanishingly small with increasing site number n. 
The relative efficiency of power transmission in the two 




FIG. 9: Time dependence of the current in in the n— th SRR 
for an array with iV = 10, A = 0.01, eo = 0.3, 7 = 0.001, and 
frequency (a) / = 0.888 GHz; (b) / = 1.78 GHz. The four 
different curves are for the current in the n = 1 SRR (black- 
solid), n = 4 SRR (red-dashed), n = 7 SRR (green-dotted), 
n = 10 SRR (blue-dot-dashed). 
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n 



FIG. 10: Power tranfer efficiency in different frequency re- 
gions, i.e., average power (in dBm) as a function of site num- 
ber n for an array £0 = 0.3, A = 0.01, 7 = 0.001, = 10, 
and / = 0.888 GHz (black diamonds); / = 1.197 GHz (red 
squares), / = 1.78 GHz (green triangles). 



the end of the extended array, and no reflection occurs. 
The total number of SRRs is now N' = N + Ni where 
Ni = 5, and the dissipation coefficient 7 is now a function 
of the site number n, 7(n), given by 

where fi = ^ with 7 — 0.001 as before. For the 

array with n— dependent dissipation coefficient, a typical 
average power density plot is shown in Fig. 11, where the 
power in each site is the average over a time interval of 
5000 Tq. We observe, at least for the first N sites of the 
array, a pattern similar with those seen in Figs. 7 and 8. 
That is, the appearance of the linear and the nonlinear 
transmission bands at about the same frequency inter- 
vals. We also observe that almost all of the transmitted 
power is absorbed from the five last sites of the array. 
Importantly, the narrow band coming from geometrical 
resonances in Figs. 7 and 8, which is around / = I GHz, 
appears in Fig. II as well, and moreover it seems that it 
is capable of transmitting power to the end of the array. 



frequencies used in the previous figure, are shown in Fig. 
10, along with the power transmission for a specific fre- 
quency (/ — 1.197 GHz) outside all linear and nonlinear 
bands. It is remarkable that the transmission efficiency 
is almost the same for frequencies in the two pass-bands, 
while for the other one the transmitted power practically 
vanish at the fourth site. 

In the density plots above, the right boundary was ac- 
tually a refiecting one, allowing the formation of station- 
ary states in the array. We have also used a totally ab- 
sorbing boundary by adding five more sites at the right 
end of the array, where the damping coefficient increases 
exponentially with increasing site number. Thus, the 
waves coming from the left are completely damped at 



V. DISCUSSION AND CONCLUDING 
REMARKS 

We have used a simple electric circuit model of in- 
ductively coupled resistive-inductive-capacitive circuits 
with voltage-dependent capacitance to investigate the 
transmission of power in a one-dimensional, end-driven, 
nonlinear magnetic metamaterial comprised of varactor- 
loaded SRRs. We have used both a reflecting and an 
absorbing boundary for the end that is not driven by the 
external field. The importance of such discrete models 
for describing the dynamic behaviour of SRR-based rnag- 
netic metamaterials has been recently appreciated [27|, 



9 




FIG. 11: Average power density plot Pn (in dBm) on the n— / 
plane for a varactor-loaded SRR array with absorbing right 
boundary (see text), for A = +0.01, eo = 0.297, 7 = 0.001, 
and N' = 15 [N = 10). 

although similar discrete models have been employed ear- 
lier for the study of localized excitations in those systems 

[Bill 

In the present work, we have focused on the trans- 
mission of power along the array, and especially that in 
the nonlinear regime, which results in the formation of 
nonlinear pass-bands for frequencies close to the reso- 
nances of a single varactor-loaded SRR that is described 
by the approximate model equation (jl3p . The approx- 
imate model holds fairly well for relatively low driving 
amplitudes, for which the dissipative diode current of the 
varactor can be neglected. It seems, thus that the prop- 
erties of the individual elements of the array predomi- 
nantly determine its transmission properties, at least in 
the case of weak magnetoinductive coupling investigated 
in the present work. The numerical results demonstrate 
that the nonlinear pass-bands transmit power with ef- 
ficiently comparable to that of the linear band. For the 
parameter sets used in this work, which are close to those 
in [HO] for a single varactor-loaded SRR, there are three 
nonlinear bands, both above and below the linear band. 



Their width apparently depends on the frequency, i.e., 
the bands at higher frequencies are significanlty wider 
than those at lower frequencies. 

Moreover, we have also observed the formation, ei- 
ther complete or partial, of another nonlinear pass-band 
at frequencies slightly above the linear band (at / ^ 
1 GHz). That band, which is rather narrow and can 
be observed both in Figs. 7 and 8 as well as in Fig. 11, 
seem to depend sensitively on the model parameters A, 
/, Eq (for a and /3 fixed), and mainly on the loss coef- 
ficient 7. Comparing Figs. 7 and 8, which have been 
obtained for different loss coefficient 7 = 0.001 and 0.01, 
respectively, we see that the power transmission in that 
band considerably decreases when the losses increase by 
an order of magnitude. That type of nonlinear band 
appears due to geometrical resonances of the standing 
waves that are formed in the array due to driving at 
one end. Geometrical resonances seem to offer an im- 
portant mechanism for the generation of nonlinear pass- 
bands, which can be employed for constructing transmis- 
sion lines with banded power transmission spectra with 
bands at desired frequency ranges. Similar conduct-free 
transmission lines could be formed by superconducting 
magnetic metamaterials, where the basic structural ele- 
ment (the SRR) is replaced b y it s direct superconducting 
analogue, the rf SQUID [13, IH. The rf SQUID, where 
the acronym stands for radio-frequency superconductive 
quantum interference device, is just a superconducting 
ring interrupted by a Josephson junction, which makes 
it an intrinsically nonlinear element. A one-dimensional 
array of rf SQUIDs coupled through their mutual induc- 
tances forms a nonlinear magnetoinductive transmission 
line [1^ whose power transmission properties is a matter 
of future work. 
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